set.seed(5168)
## responsibility attribution ##

	library(foreign)
	library(mnormt)
	library(ordinal)

  data <- read.dta('Denmark (2014).dta')

	model <- clmm(as.factor(Resp_attribution) ~ Prime_minister + Cabinet + Opposition + DK_Perc_seats + Prime_minister * DK_Perc_seats + Cabinet * DK_Perc_seats + Opposition*DK_Perc_seats + (1 | Party) + (1 | pid), data = data, HESS = TRUE, link = 'probit')
	summary(model)

	pars <- as.matrix(model$coefficients)
	vce <- vcov(model)
	vce <- vce[-c(12:13), ]      
	vce <- vce[, -c(12:13)]      
	
	simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

	summary(model)
	summary(simCoef)

	pmHi <- c()
	pmLo <- c()
	pmMed <- c()

	partHi <- c()
	partLo <- c()
	partMed <- c()

	supHi <- c()
	supLo <- c()
	supMed <- c()

	oppHi <- c()
	oppLo <- c()
	oppMed <- c()

  otherHi <- c()
  otherLo <- c()
  otherMed <- c()

  part_pred <- c()
  opp_pred <- c()
  
  marg_effect <- c()
  margeff_Hi <- c()
  margeff_Lo <- c()
  margeff_Med <- c()

for(i in 0:30){
  
  ## partner party 
  agg <- simCoef[, 5] * 0 + 
    simCoef[, 6] * 1 + 
    simCoef[, 7] * 0 + 
    simCoef[, 8] * i + 
    simCoef[, 9] * 0 + 
    simCoef[, 10] * i + 
    simCoef[, 11] * 0 
  
  pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
  pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
  pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
  pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
  pred5 <- 1 - pred1 - pred2 - pred3 - pred4
  
  part_pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  partHi[i + 1] <- quantile(pred, 0.975)
  partLo[i + 1] <- quantile(pred, 0.025)
  partMed[i + 1] <- quantile(pred, 0.5)
  
  ## opp party 
  agg <- simCoef[, 5] * 0 + 
    simCoef[, 6] * 0 + 
    simCoef[, 7] * 1 + 
    simCoef[, 8] * i + 
    simCoef[, 9] * 0 + 
    simCoef[, 10] * 0 + 
    simCoef[, 11] * i 
  
  pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
  pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
  pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
  pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
  pred5 <- 1 - pred1 - pred2 - pred3 - pred4
  
  opp_pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  oppHi[i + 1] <- quantile(pred, 0.975)
  oppLo[i + 1] <- quantile(pred, 0.025)
  oppMed[i + 1] <- quantile(pred, 0.5)
  
  #get marginal effect 
  marg_effect = part_pred - opp_pred
  margeff_Hi[i + 1] = quantile(marg_effect, 0.975)
  margeff_Lo[i + 1] = quantile(marg_effect, 0.025)
  margeff_Med[i + 1] = quantile(marg_effect, 0.5)
  
}  

#seats sequence
seats <- seq(0, 30, 1)

#Now create the figure
par(mar=c(3.5,6.2,1.5,1.5))      
plot(1, 1, type = 'n',
     xlim = c(0, 30),
     ylim = c(0, 1),
     xlab = '',
     ylab = '',
     main = '', 
     axes = FALSE)
axis(1,at=seq(0,30,10),labels=T)          
axis(2,at=seq(0,1,0.50),labels=T)
polygon(c(seats, rev(seats)), c(margeff_Hi, rev(margeff_Lo)), col = rgb(red = 0.1, 0.1, 0, 0.2), border = FALSE)
lines(seats, margeff_Med, col = rgb(0.1, 0.1, 0, 1), lwd = 3)
abline(h=0,col="red")

mtext("Change in \nresponsibility attribution", side=2, line=2.5, cex=1.2)         
mtext("Perceived seat %", side=1, line=2.2, cex=1.2)                  

dev.copy(png,'Figure 4b DK 2014 - Marginal effects of party types over seat share.png', height = 395, width = 517)  
dev.off()  
